1. Introdução

Será feita a modelagem da base AmesHousing com o pacote tidymodels e alguns dos pacotes apresentados em aula ao longo do curso de Modelos Preditivos do Programa Avançado em Data Science e Decisão do Insper

Objetivo:

Modelar uma previsão do preço de venda das casas com a maior acurácia possível da base AmesHousing.

Modelos:

  • Linear (com ou sem seleção stepwise)

  • LASSO

  • Ridge Regression

  • Bagging

  • Floresta Aleatória

Bibliotecas

library(AmesHousing)
library(tidyverse)
library(tidymodels)
library(skimr)
library(naniar)
library(GGally)
library(vip)

Dados

Será carregada a base “make_ordinal_names”, pois que nela algumas colunas já possuem fatores ordenados, o que facilitará a modelagem. Com skim, temos uma ideia do que há na base.

dados <- make_ordinal_ames()
skim(dados)
Data summary
Name dados
Number of rows 2930
Number of columns 81
_______________________
Column type frequency:
factor 46
numeric 35
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
MS_SubClass 0 1 FALSE 16 One: 1079, Two: 575, One: 287, One: 192
MS_Zoning 0 1 FALSE 7 Res: 2273, Res: 462, Flo: 139, Res: 27
Street 0 1 FALSE 2 Pav: 2918, Grv: 12
Alley 0 1 FALSE 3 No_: 2732, Gra: 120, Pav: 78
Lot_Shape 0 1 TRUE 4 Reg: 1859, Sli: 979, Mod: 76, Irr: 16
Land_Contour 0 1 TRUE 4 Lvl: 2633, HLS: 120, Bnk: 117, Low: 60
Utilities 0 1 TRUE 3 All: 2927, NoS: 2, NoS: 1, ELO: 0
Lot_Config 0 1 FALSE 5 Ins: 2140, Cor: 511, Cul: 180, FR2: 85
Land_Slope 0 1 TRUE 3 Gtl: 2789, Mod: 125, Sev: 16
Neighborhood 0 1 FALSE 28 Nor: 443, Col: 267, Old: 239, Edw: 194
Condition_1 0 1 FALSE 9 Nor: 2522, Fee: 164, Art: 92, RRA: 50
Condition_2 0 1 FALSE 8 Nor: 2900, Fee: 13, Art: 5, Pos: 4
Bldg_Type 0 1 FALSE 5 One: 2425, Twn: 233, Dup: 109, Twn: 101
House_Style 0 1 FALSE 8 One: 1481, Two: 873, One: 314, SLv: 128
Overall_Qual 0 1 TRUE 10 Ave: 825, Abo: 732, Goo: 602, Ver: 350
Overall_Cond 0 1 TRUE 9 Ave: 1654, Abo: 533, Goo: 390, Ver: 144
Roof_Style 0 1 FALSE 6 Gab: 2321, Hip: 551, Gam: 22, Fla: 20
Roof_Matl 0 1 FALSE 8 Com: 2887, Tar: 23, WdS: 9, WdS: 7
Exterior_1st 0 1 FALSE 16 Vin: 1026, Met: 450, HdB: 442, Wd : 420
Exterior_2nd 0 1 FALSE 17 Vin: 1015, Met: 447, HdB: 406, Wd : 397
Mas_Vnr_Type 0 1 FALSE 5 Non: 1775, Brk: 880, Sto: 249, Brk: 25
Exter_Qual 0 1 TRUE 4 Typ: 1799, Goo: 989, Exc: 107, Fai: 35
Exter_Cond 0 1 TRUE 5 Typ: 2549, Goo: 299, Fai: 67, Exc: 12
Foundation 0 1 FALSE 6 PCo: 1310, CBl: 1244, Brk: 311, Sla: 49
Bsmt_Qual 0 1 TRUE 6 Typ: 1283, Goo: 1219, Exc: 258, Fai: 88
Bsmt_Cond 0 1 TRUE 6 Typ: 2616, Goo: 122, Fai: 104, No_: 80
Bsmt_Exposure 0 1 TRUE 5 No: 1906, Av: 418, Gd: 284, Mn: 239
BsmtFin_Type_1 0 1 TRUE 7 GLQ: 859, Unf: 851, ALQ: 429, Rec: 288
BsmtFin_Type_2 0 1 TRUE 7 Unf: 2499, Rec: 106, LwQ: 89, No_: 81
Heating 0 1 FALSE 6 Gas: 2885, Gas: 27, Gra: 9, Wal: 6
Heating_QC 0 1 TRUE 5 Exc: 1495, Typ: 864, Goo: 476, Fai: 92
Central_Air 0 1 FALSE 2 Y: 2734, N: 196
Electrical 1 1 TRUE 5 SBr: 2682, Fus: 188, Fus: 50, Fus: 8
Kitchen_Qual 0 1 TRUE 5 Typ: 1494, Goo: 1160, Exc: 205, Fai: 70
Functional 0 1 TRUE 8 Typ: 2728, Min: 70, Min: 65, Mod: 35
Fireplace_Qu 0 1 TRUE 6 No_: 1422, Goo: 744, Typ: 600, Fai: 75
Garage_Type 0 1 FALSE 7 Att: 1731, Det: 782, Bui: 186, No_: 157
Garage_Finish 0 1 TRUE 4 Unf: 1231, RFn: 812, Fin: 728, No_: 159
Garage_Qual 0 1 TRUE 6 Typ: 2615, No_: 159, Fai: 124, Goo: 24
Garage_Cond 0 1 TRUE 6 Typ: 2665, No_: 159, Fai: 74, Goo: 15
Paved_Drive 0 1 TRUE 3 Pav: 2652, Dir: 216, Par: 62
Pool_QC 0 1 TRUE 5 No_: 2917, Goo: 4, Exc: 4, Typ: 3
Fence 0 1 TRUE 5 No_: 2358, Min: 330, Goo: 118, Goo: 112
Misc_Feature 0 1 FALSE 6 Non: 2824, She: 95, Gar: 5, Oth: 4
Sale_Type 0 1 FALSE 10 WD : 2536, New: 239, COD: 87, Con: 26
Sale_Condition 0 1 FALSE 6 Nor: 2413, Par: 245, Abn: 190, Fam: 46

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Lot_Frontage 0 1 57.65 33.50 0.00 43.00 63.00 78.00 313.00 ▇▇▁▁▁
Lot_Area 0 1 10147.92 7880.02 1300.00 7440.25 9436.50 11555.25 215245.00 ▇▁▁▁▁
Year_Built 0 1 1971.36 30.25 1872.00 1954.00 1973.00 2001.00 2010.00 ▁▂▃▆▇
Year_Remod_Add 0 1 1984.27 20.86 1950.00 1965.00 1993.00 2004.00 2010.00 ▅▂▂▃▇
Mas_Vnr_Area 0 1 101.10 178.63 0.00 0.00 0.00 162.75 1600.00 ▇▁▁▁▁
BsmtFin_SF_1 0 1 4.18 2.23 0.00 3.00 3.00 7.00 7.00 ▃▂▇▁▇
BsmtFin_SF_2 0 1 49.71 169.14 0.00 0.00 0.00 0.00 1526.00 ▇▁▁▁▁
Bsmt_Unf_SF 0 1 559.07 439.54 0.00 219.00 465.50 801.75 2336.00 ▇▅▂▁▁
Total_Bsmt_SF 0 1 1051.26 440.97 0.00 793.00 990.00 1301.50 6110.00 ▇▃▁▁▁
First_Flr_SF 0 1 1159.56 391.89 334.00 876.25 1084.00 1384.00 5095.00 ▇▃▁▁▁
Second_Flr_SF 0 1 335.46 428.40 0.00 0.00 0.00 703.75 2065.00 ▇▃▂▁▁
Low_Qual_Fin_SF 0 1 4.68 46.31 0.00 0.00 0.00 0.00 1064.00 ▇▁▁▁▁
Gr_Liv_Area 0 1 1499.69 505.51 334.00 1126.00 1442.00 1742.75 5642.00 ▇▇▁▁▁
Bsmt_Full_Bath 0 1 0.43 0.52 0.00 0.00 0.00 1.00 3.00 ▇▆▁▁▁
Bsmt_Half_Bath 0 1 0.06 0.25 0.00 0.00 0.00 0.00 2.00 ▇▁▁▁▁
Full_Bath 0 1 1.57 0.55 0.00 1.00 2.00 2.00 4.00 ▁▇▇▁▁
Half_Bath 0 1 0.38 0.50 0.00 0.00 0.00 1.00 2.00 ▇▁▅▁▁
Bedroom_AbvGr 0 1 2.85 0.83 0.00 2.00 3.00 3.00 8.00 ▁▇▂▁▁
Kitchen_AbvGr 0 1 1.04 0.21 0.00 1.00 1.00 1.00 3.00 ▁▇▁▁▁
TotRms_AbvGrd 0 1 6.44 1.57 2.00 5.00 6.00 7.00 15.00 ▁▇▂▁▁
Fireplaces 0 1 0.60 0.65 0.00 0.00 1.00 1.00 4.00 ▇▇▁▁▁
Garage_Cars 0 1 1.77 0.76 0.00 1.00 2.00 2.00 5.00 ▅▇▂▁▁
Garage_Area 0 1 472.66 215.19 0.00 320.00 480.00 576.00 1488.00 ▃▇▃▁▁
Wood_Deck_SF 0 1 93.75 126.36 0.00 0.00 0.00 168.00 1424.00 ▇▁▁▁▁
Open_Porch_SF 0 1 47.53 67.48 0.00 0.00 27.00 70.00 742.00 ▇▁▁▁▁
Enclosed_Porch 0 1 23.01 64.14 0.00 0.00 0.00 0.00 1012.00 ▇▁▁▁▁
Three_season_porch 0 1 2.59 25.14 0.00 0.00 0.00 0.00 508.00 ▇▁▁▁▁
Screen_Porch 0 1 16.00 56.09 0.00 0.00 0.00 0.00 576.00 ▇▁▁▁▁
Pool_Area 0 1 2.24 35.60 0.00 0.00 0.00 0.00 800.00 ▇▁▁▁▁
Misc_Val 0 1 50.64 566.34 0.00 0.00 0.00 0.00 17000.00 ▇▁▁▁▁
Mo_Sold 0 1 6.22 2.71 1.00 4.00 6.00 8.00 12.00 ▅▆▇▃▃
Year_Sold 0 1 2007.79 1.32 2006.00 2007.00 2008.00 2009.00 2010.00 ▇▇▇▇▃
Sale_Price 0 1 180796.06 79886.69 12789.00 129500.00 160000.00 213500.00 755000.00 ▇▇▁▁▁
Longitude 0 1 -93.64 0.03 -93.69 -93.66 -93.64 -93.62 -93.58 ▅▅▇▆▁
Latitude 0 1 42.03 0.02 41.99 42.02 42.03 42.05 42.06 ▂▂▇▇▇

Aparentemente existe um valor faltando na coluna Electrical. Vamos carregar a base make_ames e verificar se esse valor está faltando também, novamnte com skim.

dados2 <- make_ames()
skim(dados2$Electrical)
Data summary
Name dados2$Electrical
Number of rows 2930
Number of columns 1
_______________________
Column type frequency:
factor 1
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
data 0 1 FALSE 6 SBr: 2682, Fus: 188, Fus: 50, Fus: 8

Pelo visto não há nenhum valor faltando nessa coluna. Vamos olhar mais de perto. Primeiro encontramos a linha com valor faltando na base original, e depois procuramos o valor dessa linha na nova base.

row_na <- dados %>% 
  rowid_to_column() %>%
  filter(is.na(dados$Electrical)) %>% 
  select(rowid)

(dados$Electrical[row_na$rowid])
## [1] <NA>
## Levels: Mix < FuseP < FuseF < FuseA < SBrkr
(dados2$Electrical[row_na$rowid])
## [1] Unknown
## Levels: FuseA FuseF FuseP Mix SBrkr Unknown

O valor é “Unknown”, assim como na base ordenada. Como em ambas não há a informação de “Electrical” em apenas uma linha, vamos utilizar dropna() na base ordenada e seguir com a análise.

dados <- dados %>% 
  drop_na()

skim(dados$Electrical)
Data summary
Name dados$Electrical
Number of rows 2929
Number of columns 1
_______________________
Column type frequency:
factor 1
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
data 0 1 TRUE 5 SBr: 2682, Fus: 188, Fus: 50, Fus: 8

2. Análise Exploratória

Para entender melhor o compartamento dos dados, será feita uma análise exploratória da base buscando relações com a variável predita “Sale_Price”.

2.1. Distribuição dos preços

Iniciamos a análise entendendo o perfil dessa variável:

  • Por região
dados %>%
  ggplot(aes(Longitude, Latitude, color = Sale_Price)) +
  geom_point(size = 0.5, alpha = 0.7) +
  scale_color_gradient(low = 'skyblue', high = 'darkred', labels = scales::label_number_si())

Verifica-se uma concentração de casas mais caras principalmente no norte da região analisada. Provavelmente ‘Neighborhood’ (bairro), assim como Latitude e Longitude, devem ter uma correlação alta com o preço.

  • Bairro
dados %>%
  ggplot(aes(Neighborhood, Sale_Price, fill = Neighborhood)) +
  geom_boxplot(show.legend = FALSE)+
  scale_y_continuous(labels = scales::label_number_si())+
  coord_flip()

  • Histograma
dados %>% 
    ggplot(aes(x=Sale_Price)) +
    geom_histogram(fill="skyblue", binwidth = 10000) +
    scale_x_continuous(breaks= seq(0, 800000, by=100000), labels = scales::label_number_si())

(summary(dados$Sale_Price))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12789  129500  160000  180801  213500  755000

No histograma observa-se uma concentração mais perto do começo da amostra. Isso é esperado, pois menos pessoas conseguem comprar casas mais caras. Importante notar também que a média e mediana são próximas.

2.2. Correlação entre variáveis numéricas

Nesse item o objetivo é ter uma ideia de quais variáveis estão mais correlacionadas com “Sale_Price” para ter uma ideia do preditores e para nortear o restante da análise exploratória.

 #variáveis numéricas
num_vars <- which(sapply(dados, is.numeric))

#correlação entre todas as variáveis
corr_all <- cor(dados[,num_vars], use="pairwise.complete.obs") 

#correlação com Sale_Price em ordem decrescente
(corr_sorted <- as.matrix(sort(corr_all[,'Sale_Price'], decreasing = TRUE))) 
##                           [,1]
## Sale_Price          1.00000000
## Gr_Liv_Area         0.70677666
## Garage_Cars         0.64759257
## Garage_Area         0.64013460
## Total_Bsmt_SF       0.63269326
## First_Flr_SF        0.62173389
## Year_Built          0.55861903
## Full_Bath           0.54570831
## Year_Remod_Add      0.53314636
## Mas_Vnr_Area        0.50219365
## TotRms_AbvGrd       0.49550750
## Fireplaces          0.47457710
## Wood_Deck_SF        0.32714767
## Open_Porch_SF       0.31293846
## Latitude            0.29101035
## Half_Bath           0.28520178
## Bsmt_Full_Bath      0.27580905
## Second_Flr_SF       0.26943829
## Lot_Area            0.26654763
## Lot_Frontage        0.20190876
## Bsmt_Unf_SF         0.18329078
## Bedroom_AbvGr       0.14392488
## Screen_Porch        0.11213709
## Pool_Area           0.06840003
## Mo_Sold             0.03523475
## Three_season_porch  0.03221900
## BsmtFin_SF_2        0.00600098
## Misc_Val           -0.01569664
## Year_Sold          -0.03056032
## Bsmt_Half_Bath     -0.03583132
## Low_Qual_Fin_SF    -0.03766575
## Kitchen_AbvGr      -0.11982695
## Enclosed_Porch     -0.12881128
## BsmtFin_SF_1       -0.13487107
## Longitude          -0.25141580
#correlação alta com Sale_Price
corr_high <- names(which(apply(corr_sorted, 1, function(x) abs(x)>0.5))) 

dados[,corr_high] %>% 
  ggpairs()

2.3. Ano de construção, Área total e Material de construção

dados %>%
  ggplot(aes(Year_Built, Sale_Price, size = Lot_Area/1000,  color = Foundation)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

O mais marcante nesse gráfico é a mudança do material utilizado na estrutura das casas ao longo dos anos. A relação entre ano de construção e preço é também confirmada, visto que ‘Year_Built’ é umas das variáveis numéricas com correlação maior do que 0.5.

2.4. Garagem

No topo da lista de variáveis correlacionadas está a a área da garagem e quantidade de carros. Vamos observar o comportamento destas e outras variáveis relacionadas.

dados %>%
  ggplot(aes(Garage_Area, Sale_Price, size = Garage_Cars, color = Garage_Type)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

dados %>%
  ggplot(aes(Garage_Qual, Sale_Price, color = Garage_Finish)) +
  geom_boxplot()+
  scale_y_continuous(labels = scales::label_number_si())

Além das variáveis altamente correlacioanadas, há uma concentração do tipo de garagem em difenrentes níveis de preço. O tipo de acabamento, no entanto, parece não ter nenhuma relação significativa.

2.5. Área total, Lareiras e Condição de Venda

dados %>%
  ggplot(aes(Gr_Liv_Area, Sale_Price, size = Fireplaces, color = Sale_Condition)) +
  geom_point(alpha = 0.5)+
  scale_y_continuous(labels = scales::label_number_si())

Aqui apenas é confirmada a correlação da quantidade de lareiras e área comum. A condição de venda “Normal” predomina na amostra.

2.5 Qualidade de venda

dados %>%
  ggplot(aes(Overall_Qual, Sale_Price, fill = Street)) +
  geom_boxplot()+
  scale_y_continuous(labels = scales::label_number_si())

Apesar de estar na base como uma variável categórica, “Overall Quality” parace ter uma correlação alta com Sale_Price.

3. Pré processamento

Antes de alimentar os dados aos modelos, é necessário separar a base em treino e teste, e garantir que está em um formato interpretável para os algoritmos utilizados. O parâmetro “strata” divide a base proporcionalmente conforme a variável indicada.

3.1 Separação da base em treino e teste (treino/teste/total):

set.seed(123)
(ames_split <- initial_split(dados, prop = 0.8, strata = 'Sale_Price'))
## <2345/584/2929>
ames_train <- training(ames_split)
ames_test <- testing(ames_split)

3.2. Receita

A interface do tidymodels permite a criação de uma receita com o pacote “recipes” para processar os dados antes do modelo. Isso facilita também o pré-processamento da base de teste ao final do relatório, pois aplica exatamente as mesmas modificações. Algumas variáveis já estão na base como fator, portanto não é necessário ordená-las. No entanto, é necessário incluir na receita um passo para converter a ordem em números: step_ordinalscore

A seguir as variáveis ordenadas:

ord_vars <- vapply(dados, is.ordered, logical(1))

(ordered <-names(ord_vars)[ord_vars])
##  [1] "Lot_Shape"      "Land_Contour"   "Utilities"      "Land_Slope"    
##  [5] "Overall_Qual"   "Overall_Cond"   "Exter_Qual"     "Exter_Cond"    
##  [9] "Bsmt_Qual"      "Bsmt_Cond"      "Bsmt_Exposure"  "BsmtFin_Type_1"
## [13] "BsmtFin_Type_2" "Heating_QC"     "Electrical"     "Kitchen_Qual"  
## [17] "Functional"     "Fireplace_Qu"   "Garage_Finish"  "Garage_Qual"   
## [21] "Garage_Cond"    "Paved_Drive"    "Pool_QC"        "Fence"

Como algumas colunas possuem muitas categorias, as que possuem mais de 10 serão reduzidas com “step_other” (10 é o número de classificações atribuídas nas colunas que avaliam a qualidade). Após isso, cria-se variáveis dummy para as demais categorias (step_dummy), remove-se as variáveis com apenas um valor (step_vz) e normaliza-se a base (step_normalize).

(ames_rec <- recipe(Sale_Price ~ ., data = ames_train) %>%
  #step_log(Sale_Price)) %>% 
  step_other(MS_SubClass, Neighborhood, Exterior_1st, Exterior_2nd, threshold = 0.02) %>% 
  step_ordinalscore(ordered) %>%
  step_dummy(all_nominal(), -all_outcomes()) %>%
  step_zv(all_predictors()) %>% 
  step_normalize(all_numeric()) %>% 
  prep())
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         80
## 
## Training data contained 2345 data points and no missing data.
## 
## Operations:
## 
## Collapsing factor levels for MS_SubClass, Neighborhood, ... [trained]
## Scoring for Lot_Shape, Land_Contour, Utilities, ... [trained]
## Dummy variables from MS_SubClass, MS_Zoning, Street, Alley, ... [trained]
## Zero variance filter removed Roof_Matl_Membran [trained]
## Centering and scaling for Lot_Frontage, Lot_Area, Lot_Shape, ... [trained]

A receita pronta é então aplicada nas bases de treino e teste com a função “bake”. Para a base teste, pode-se utilizar juice, que é um caso especial de bake.

train_baked  <- juice(ames_rec) 

test_baked <- bake(ames_rec, new_data = ames_test)

A seguir uma ideia de como ficou a base processada:

skim(train_baked)
Data summary
Name train_baked
Number of rows 2345
Number of columns 190
_______________________
Column type frequency:
numeric 190
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Lot_Frontage 0 1 0 1 -1.75 -0.42 0.15 0.64 7.73 ▇▇▁▁▁
Lot_Area 0 1 0 1 -1.10 -0.33 -0.08 0.19 25.71 ▇▁▁▁▁
Lot_Shape 0 1 0 1 -4.61 -1.07 0.70 0.70 0.70 ▁▁▁▅▇
Land_Contour 0 1 0 1 -4.67 0.31 0.31 0.31 0.31 ▁▁▁▁▇
Utilities 0 1 0 1 -34.22 0.03 0.03 0.03 0.03 ▁▁▁▁▇
Land_Slope 0 1 0 1 -7.87 0.22 0.22 0.22 0.22 ▁▁▁▁▇
Overall_Qual 0 1 0 1 -3.62 -0.77 -0.06 0.65 2.79 ▁▂▇▅▁
Overall_Cond 0 1 0 1 -4.18 -0.52 -0.52 0.40 3.15 ▁▁▇▅▁
Year_Built 0 1 0 1 -3.25 -0.59 0.06 0.95 1.28 ▁▂▃▆▇
Year_Remod_Add 0 1 0 1 -1.64 -0.92 0.43 0.91 1.24 ▅▂▂▃▇
Mas_Vnr_Area 0 1 0 1 -0.57 -0.57 -0.57 0.36 6.84 ▇▁▁▁▁
Exter_Qual 0 1 0 1 -2.40 -0.68 -0.68 1.05 2.77 ▁▇▁▅▁
Exter_Cond 0 1 0 1 -5.72 -0.23 -0.23 -0.23 5.26 ▁▁▇▁▁
Bsmt_Qual 0 1 0 1 -3.87 -0.53 0.58 0.58 1.69 ▁▁▇▇▂
Bsmt_Cond 0 1 0 1 -5.14 0.12 0.12 0.12 3.63 ▁▁▇▁▁
Bsmt_Exposure 0 1 0 1 -1.52 -0.58 -0.58 0.35 2.23 ▁▇▁▂▁
BsmtFin_Type_1 0 1 0 1 -1.69 -1.22 0.20 1.15 1.15 ▆▁▂▂▇
BsmtFin_SF_1 0 1 0 1 -1.41 -0.52 -0.52 1.27 1.27 ▅▆▁▁▇
BsmtFin_Type_2 0 1 0 1 -1.33 -0.30 -0.30 -0.30 4.86 ▇▁▁▁▁
BsmtFin_SF_2 0 1 0 1 -0.30 -0.30 -0.30 -0.30 8.58 ▇▁▁▁▁
Bsmt_Unf_SF 0 1 0 1 -1.27 -0.77 -0.21 0.55 4.07 ▇▅▂▁▁
Total_Bsmt_SF 0 1 0 1 -2.46 -0.61 -0.14 0.58 9.47 ▇▇▁▁▁
Heating_QC 0 1 0 1 -3.29 -1.20 0.89 0.89 0.89 ▁▁▅▂▇
Electrical 0 1 0 1 -9.89 0.28 0.28 0.28 0.28 ▁▁▁▁▇
First_Flr_SF 0 1 0 1 -2.13 -0.73 -0.19 0.59 10.22 ▇▃▁▁▁
Second_Flr_SF 0 1 0 1 -0.79 -0.79 -0.79 0.85 4.04 ▇▃▂▁▁
Low_Qual_Fin_SF 0 1 0 1 -0.10 -0.10 -0.10 -0.10 22.30 ▇▁▁▁▁
Gr_Liv_Area 0 1 0 1 -2.35 -0.75 -0.12 0.49 7.26 ▅▇▁▁▁
Bsmt_Full_Bath 0 1 0 1 -0.82 -0.82 -0.82 1.08 4.88 ▇▆▁▁▁
Bsmt_Half_Bath 0 1 0 1 -0.25 -0.25 -0.25 -0.25 7.84 ▇▁▁▁▁
Full_Bath 0 1 0 1 -2.81 -1.01 0.78 0.78 4.38 ▁▇▇▁▁
Half_Bath 0 1 0 1 -0.76 -0.76 -0.76 1.22 3.19 ▇▁▅▁▁
Bedroom_AbvGr 0 1 0 1 -3.41 -1.02 0.18 0.18 6.16 ▁▇▂▁▁
Kitchen_AbvGr 0 1 0 1 -4.92 -0.20 -0.20 -0.20 9.22 ▁▇▁▁▁
Kitchen_Qual 0 1 0 1 -2.28 -0.76 -0.76 0.76 2.28 ▁▇▁▆▁
TotRms_AbvGrd 0 1 0 1 -2.80 -0.91 -0.28 0.35 5.39 ▁▇▂▁▁
Functional 0 1 0 1 -10.53 0.24 0.24 0.24 0.24 ▁▁▁▁▇
Fireplaces 0 1 0 1 -0.93 -0.93 0.61 0.61 5.23 ▇▇▁▁▁
Fireplace_Qu 0 1 0 1 -0.99 -0.99 0.11 1.22 1.77 ▇▁▃▅▁
Garage_Finish 0 1 0 1 -1.91 -0.79 0.32 0.32 1.44 ▁▇▁▆▅
Garage_Cars 0 1 0 1 -2.31 -1.00 0.31 0.31 4.23 ▁▃▇▂▁
Garage_Area 0 1 0 1 -2.20 -0.70 0.04 0.49 4.75 ▃▇▂▁▁
Garage_Qual 0 1 0 1 -3.87 0.28 0.28 0.28 3.05 ▁▁▁▇▁
Garage_Cond 0 1 0 1 -3.90 0.27 0.27 0.27 3.05 ▁▁▁▇▁
Paved_Drive 0 1 0 1 -3.50 0.31 0.31 0.31 0.31 ▁▁▁▁▇
Wood_Deck_SF 0 1 0 1 -0.75 -0.75 -0.75 0.61 6.26 ▇▂▁▁▁
Open_Porch_SF 0 1 0 1 -0.70 -0.70 -0.32 0.33 10.20 ▇▁▁▁▁
Enclosed_Porch 0 1 0 1 -0.38 -0.38 -0.38 -0.38 8.72 ▇▁▁▁▁
Three_season_porch 0 1 0 1 -0.10 -0.10 -0.10 -0.10 19.11 ▇▁▁▁▁
Screen_Porch 0 1 0 1 -0.29 -0.29 -0.29 -0.29 9.92 ▇▁▁▁▁
Pool_Area 0 1 0 1 -0.06 -0.06 -0.06 -0.06 21.31 ▇▁▁▁▁
Pool_QC 0 1 0 1 -0.07 -0.07 -0.07 -0.07 18.33 ▇▁▁▁▁
Fence 0 1 0 1 -0.48 -0.48 -0.48 -0.48 2.76 ▇▁▁▁▁
Misc_Val 0 1 0 1 -0.09 -0.09 -0.09 -0.09 27.08 ▇▁▁▁▁
Mo_Sold 0 1 0 1 -1.92 -0.82 -0.08 0.65 2.13 ▃▆▇▃▃
Year_Sold 0 1 0 1 -1.36 -0.59 0.17 0.93 1.69 ▇▇▇▇▃
Longitude 0 1 0 1 -1.94 -0.67 0.04 0.80 2.56 ▅▅▇▆▁
Latitude 0 1 0 1 -2.63 -0.68 0.01 0.83 1.58 ▂▂▇▇▇
Sale_Price 0 1 0 1 -2.09 -0.64 -0.26 0.40 7.14 ▇▇▁▁▁
MS_SubClass_One_Story_1945_and_Older 0 1 0 1 -0.22 -0.22 -0.22 -0.22 4.55 ▇▁▁▁▁
MS_SubClass_One_and_Half_Story_Finished_All_Ages 0 1 0 1 -0.34 -0.34 -0.34 -0.34 2.97 ▇▁▁▁▁
MS_SubClass_Two_Story_1946_and_Newer 0 1 0 1 -0.50 -0.50 -0.50 -0.50 2.01 ▇▁▁▁▂
MS_SubClass_Two_Story_1945_and_Older 0 1 0 1 -0.22 -0.22 -0.22 -0.22 4.57 ▇▁▁▁▁
MS_SubClass_Split_or_Multilevel 0 1 0 1 -0.20 -0.20 -0.20 -0.20 5.09 ▇▁▁▁▁
MS_SubClass_Duplex_All_Styles_and_Ages 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.19 ▇▁▁▁▁
MS_SubClass_One_Story_PUD_1946_and_Newer 0 1 0 1 -0.27 -0.27 -0.27 -0.27 3.75 ▇▁▁▁▁
MS_SubClass_Two_Story_PUD_1946_and_Newer 0 1 0 1 -0.21 -0.21 -0.21 -0.21 4.66 ▇▁▁▁▁
MS_SubClass_Two_Family_conversion_All_Styles_and_Ages 0 1 0 1 -0.15 -0.15 -0.15 -0.15 6.84 ▇▁▁▁▁
MS_SubClass_other 0 1 0 1 -0.20 -0.20 -0.20 -0.20 4.95 ▇▁▁▁▁
MS_Zoning_Residential_High_Density 0 1 0 1 -0.10 -0.10 -0.10 -0.10 10.27 ▇▁▁▁▁
MS_Zoning_Residential_Low_Density 0 1 0 1 -1.85 0.54 0.54 0.54 0.54 ▂▁▁▁▇
MS_Zoning_Residential_Medium_Density 0 1 0 1 -0.43 -0.43 -0.43 -0.43 2.31 ▇▁▁▁▂
MS_Zoning_A_agr 0 1 0 1 -0.03 -0.03 -0.03 -0.03 34.22 ▇▁▁▁▁
MS_Zoning_C_all 0 1 0 1 -0.10 -0.10 -0.10 -0.10 10.27 ▇▁▁▁▁
MS_Zoning_I_all 0 1 0 1 -0.03 -0.03 -0.03 -0.03 34.22 ▇▁▁▁▁
Street_Pave 0 1 0 1 -15.28 0.07 0.07 0.07 0.07 ▁▁▁▁▇
Alley_No_Alley_Access 0 1 0 1 -3.65 0.27 0.27 0.27 0.27 ▁▁▁▁▇
Alley_Paved 0 1 0 1 -0.16 -0.16 -0.16 -0.16 6.12 ▇▁▁▁▁
Lot_Config_CulDSac 0 1 0 1 -0.26 -0.26 -0.26 -0.26 3.89 ▇▁▁▁▁
Lot_Config_FR2 0 1 0 1 -0.17 -0.17 -0.17 -0.17 5.83 ▇▁▁▁▁
Lot_Config_FR3 0 1 0 1 -0.07 -0.07 -0.07 -0.07 15.28 ▇▁▁▁▁
Lot_Config_Inside 0 1 0 1 -1.65 -1.65 0.61 0.61 0.61 ▃▁▁▁▇
Neighborhood_College_Creek 0 1 0 1 -0.32 -0.32 -0.32 -0.32 3.15 ▇▁▁▁▁
Neighborhood_Old_Town 0 1 0 1 -0.30 -0.30 -0.30 -0.30 3.37 ▇▁▁▁▁
Neighborhood_Edwards 0 1 0 1 -0.27 -0.27 -0.27 -0.27 3.71 ▇▁▁▁▁
Neighborhood_Somerset 0 1 0 1 -0.26 -0.26 -0.26 -0.26 3.82 ▇▁▁▁▁
Neighborhood_Northridge_Heights 0 1 0 1 -0.24 -0.24 -0.24 -0.24 4.11 ▇▁▁▁▁
Neighborhood_Gilbert 0 1 0 1 -0.24 -0.24 -0.24 -0.24 4.18 ▇▁▁▁▁
Neighborhood_Sawyer 0 1 0 1 -0.23 -0.23 -0.23 -0.23 4.38 ▇▁▁▁▁
Neighborhood_Northwest_Ames 0 1 0 1 -0.22 -0.22 -0.22 -0.22 4.57 ▇▁▁▁▁
Neighborhood_Sawyer_West 0 1 0 1 -0.22 -0.22 -0.22 -0.22 4.59 ▇▁▁▁▁
Neighborhood_Mitchell 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.19 ▇▁▁▁▁
Neighborhood_Brookside 0 1 0 1 -0.20 -0.20 -0.20 -0.20 4.89 ▇▁▁▁▁
Neighborhood_Crawford 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.22 ▇▁▁▁▁
Neighborhood_Iowa_DOT_and_Rail_Road 0 1 0 1 -0.18 -0.18 -0.18 -0.18 5.54 ▇▁▁▁▁
Neighborhood_Timberland 0 1 0 1 -0.16 -0.16 -0.16 -0.16 6.39 ▇▁▁▁▁
Neighborhood_Northridge 0 1 0 1 -0.15 -0.15 -0.15 -0.15 6.45 ▇▁▁▁▁
Neighborhood_other 0 1 0 1 -0.34 -0.34 -0.34 -0.34 2.90 ▇▁▁▁▁
Condition_1_Feedr 0 1 0 1 -0.25 -0.25 -0.25 -0.25 4.03 ▇▁▁▁▁
Condition_1_Norm 0 1 0 1 -2.49 0.40 0.40 0.40 0.40 ▁▁▁▁▇
Condition_1_PosA 0 1 0 1 -0.08 -0.08 -0.08 -0.08 12.46 ▇▁▁▁▁
Condition_1_PosN 0 1 0 1 -0.11 -0.11 -0.11 -0.11 8.78 ▇▁▁▁▁
Condition_1_RRAe 0 1 0 1 -0.09 -0.09 -0.09 -0.09 11.06 ▇▁▁▁▁
Condition_1_RRAn 0 1 0 1 -0.14 -0.14 -0.14 -0.14 7.40 ▇▁▁▁▁
Condition_1_RRNe 0 1 0 1 -0.05 -0.05 -0.05 -0.05 21.63 ▇▁▁▁▁
Condition_1_RRNn 0 1 0 1 -0.05 -0.05 -0.05 -0.05 19.74 ▇▁▁▁▁
Condition_2_Feedr 0 1 0 1 -0.07 -0.07 -0.07 -0.07 15.28 ▇▁▁▁▁
Condition_2_Norm 0 1 0 1 -10.78 0.09 0.09 0.09 0.09 ▁▁▁▁▇
Condition_2_PosA 0 1 0 1 -0.04 -0.04 -0.04 -0.04 27.93 ▇▁▁▁▁
Condition_2_PosN 0 1 0 1 -0.03 -0.03 -0.03 -0.03 34.22 ▇▁▁▁▁
Condition_2_RRAe 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Condition_2_RRAn 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Condition_2_RRNn 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Bldg_Type_TwoFmCon 0 1 0 1 -0.15 -0.15 -0.15 -0.15 6.71 ▇▁▁▁▁
Bldg_Type_Duplex 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.19 ▇▁▁▁▁
Bldg_Type_Twnhs 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.29 ▇▁▁▁▁
Bldg_Type_TwnhsE 0 1 0 1 -0.30 -0.30 -0.30 -0.30 3.37 ▇▁▁▁▁
House_Style_One_and_Half_Unf 0 1 0 1 -0.09 -0.09 -0.09 -0.09 11.70 ▇▁▁▁▁
House_Style_One_Story 0 1 0 1 -1.00 -1.00 1.00 1.00 1.00 ▇▁▁▁▇
House_Style_SFoyer 0 1 0 1 -0.17 -0.17 -0.17 -0.17 5.83 ▇▁▁▁▁
House_Style_SLvl 0 1 0 1 -0.20 -0.20 -0.20 -0.20 4.89 ▇▁▁▁▁
House_Style_Two_and_Half_Fin 0 1 0 1 -0.05 -0.05 -0.05 -0.05 19.74 ▇▁▁▁▁
House_Style_Two_and_Half_Unf 0 1 0 1 -0.09 -0.09 -0.09 -0.09 10.78 ▇▁▁▁▁
House_Style_Two_Story 0 1 0 1 -0.65 -0.65 -0.65 1.53 1.53 ▇▁▁▁▃
Roof_Style_Gable 0 1 0 1 -1.95 0.51 0.51 0.51 0.51 ▂▁▁▁▇
Roof_Style_Gambrel 0 1 0 1 -0.09 -0.09 -0.09 -0.09 11.70 ▇▁▁▁▁
Roof_Style_Hip 0 1 0 1 -0.48 -0.48 -0.48 -0.48 2.08 ▇▁▁▁▂
Roof_Style_Mansard 0 1 0 1 -0.06 -0.06 -0.06 -0.06 16.11 ▇▁▁▁▁
Roof_Style_Shed 0 1 0 1 -0.05 -0.05 -0.05 -0.05 21.63 ▇▁▁▁▁
Roof_Matl_CompShg 0 1 0 1 -8.78 0.11 0.11 0.11 0.11 ▁▁▁▁▇
Roof_Matl_Metal 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Roof_Matl_Roll 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Roof_Matl_Tar.Grv 0 1 0 1 -0.09 -0.09 -0.09 -0.09 10.78 ▇▁▁▁▁
Roof_Matl_WdShake 0 1 0 1 -0.04 -0.04 -0.04 -0.04 24.19 ▇▁▁▁▁
Roof_Matl_WdShngl 0 1 0 1 -0.04 -0.04 -0.04 -0.04 24.19 ▇▁▁▁▁
Exterior_1st_CemntBd 0 1 0 1 -0.22 -0.22 -0.22 -0.22 4.62 ▇▁▁▁▁
Exterior_1st_HdBoard 0 1 0 1 -0.41 -0.41 -0.41 -0.41 2.44 ▇▁▁▁▂
Exterior_1st_MetalSd 0 1 0 1 -0.43 -0.43 -0.43 -0.43 2.34 ▇▁▁▁▂
Exterior_1st_Plywood 0 1 0 1 -0.29 -0.29 -0.29 -0.29 3.48 ▇▁▁▁▁
Exterior_1st_VinylSd 0 1 0 1 -0.74 -0.74 -0.74 1.36 1.36 ▇▁▁▁▅
Exterior_1st_Wd.Sdng 0 1 0 1 -0.41 -0.41 -0.41 -0.41 2.42 ▇▁▁▁▂
Exterior_1st_WdShing 0 1 0 1 -0.14 -0.14 -0.14 -0.14 6.92 ▇▁▁▁▁
Exterior_1st_other 0 1 0 1 -0.19 -0.19 -0.19 -0.19 5.25 ▇▁▁▁▁
Exterior_2nd_HdBoard 0 1 0 1 -0.39 -0.39 -0.39 -0.39 2.59 ▇▁▁▁▁
Exterior_2nd_MetalSd 0 1 0 1 -0.43 -0.43 -0.43 -0.43 2.34 ▇▁▁▁▂
Exterior_2nd_Plywood 0 1 0 1 -0.32 -0.32 -0.32 -0.32 3.10 ▇▁▁▁▁
Exterior_2nd_VinylSd 0 1 0 1 -0.73 -0.73 -0.73 1.37 1.37 ▇▁▁▁▅
Exterior_2nd_Wd.Sdng 0 1 0 1 -0.40 -0.40 -0.40 -0.40 2.49 ▇▁▁▁▁
Exterior_2nd_Wd.Shng 0 1 0 1 -0.17 -0.17 -0.17 -0.17 5.87 ▇▁▁▁▁
Exterior_2nd_other 0 1 0 1 -0.26 -0.26 -0.26 -0.26 3.89 ▇▁▁▁▁
Mas_Vnr_Type_BrkFace 0 1 0 1 -0.64 -0.64 -0.64 1.55 1.55 ▇▁▁▁▃
Mas_Vnr_Type_CBlock 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Mas_Vnr_Type_None 0 1 0 1 -1.25 -1.25 0.80 0.80 0.80 ▅▁▁▁▇
Mas_Vnr_Type_Stone 0 1 0 1 -0.31 -0.31 -0.31 -0.31 3.25 ▇▁▁▁▁
Foundation_CBlock 0 1 0 1 -0.86 -0.86 -0.86 1.17 1.17 ▇▁▁▁▆
Foundation_PConc 0 1 0 1 -0.90 -0.90 -0.90 1.11 1.11 ▇▁▁▁▆
Foundation_Slab 0 1 0 1 -0.12 -0.12 -0.12 -0.12 8.12 ▇▁▁▁▁
Foundation_Stone 0 1 0 1 -0.05 -0.05 -0.05 -0.05 18.27 ▇▁▁▁▁
Foundation_Wood 0 1 0 1 -0.04 -0.04 -0.04 -0.04 24.19 ▇▁▁▁▁
Heating_GasA 0 1 0 1 -7.69 0.13 0.13 0.13 0.13 ▁▁▁▁▇
Heating_GasW 0 1 0 1 -0.10 -0.10 -0.10 -0.10 9.63 ▇▁▁▁▁
Heating_Grav 0 1 0 1 -0.05 -0.05 -0.05 -0.05 18.27 ▇▁▁▁▁
Heating_OthW 0 1 0 1 -0.03 -0.03 -0.03 -0.03 34.22 ▇▁▁▁▁
Heating_Wall 0 1 0 1 -0.04 -0.04 -0.04 -0.04 24.19 ▇▁▁▁▁
Central_Air_Y 0 1 0 1 -3.62 0.28 0.28 0.28 0.28 ▁▁▁▁▇
Garage_Type_Basment 0 1 0 1 -0.11 -0.11 -0.11 -0.11 8.93 ▇▁▁▁▁
Garage_Type_BuiltIn 0 1 0 1 -0.26 -0.26 -0.26 -0.26 3.78 ▇▁▁▁▁
Garage_Type_CarPort 0 1 0 1 -0.07 -0.07 -0.07 -0.07 14.56 ▇▁▁▁▁
Garage_Type_Detchd 0 1 0 1 -0.60 -0.60 -0.60 1.67 1.67 ▇▁▁▁▃
Garage_Type_More_Than_Two_Types 0 1 0 1 -0.09 -0.09 -0.09 -0.09 10.78 ▇▁▁▁▁
Garage_Type_No_Garage 0 1 0 1 -0.24 -0.24 -0.24 -0.24 4.14 ▇▁▁▁▁
Misc_Feature_Gar2 0 1 0 1 -0.05 -0.05 -0.05 -0.05 21.63 ▇▁▁▁▁
Misc_Feature_None 0 1 0 1 -5.09 0.20 0.20 0.20 0.20 ▁▁▁▁▇
Misc_Feature_Othr 0 1 0 1 -0.04 -0.04 -0.04 -0.04 27.93 ▇▁▁▁▁
Misc_Feature_Shed 0 1 0 1 -0.18 -0.18 -0.18 -0.18 5.43 ▇▁▁▁▁
Misc_Feature_TenC 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Sale_Type_Con 0 1 0 1 -0.04 -0.04 -0.04 -0.04 24.19 ▇▁▁▁▁
Sale_Type_ConLD 0 1 0 1 -0.10 -0.10 -0.10 -0.10 10.27 ▇▁▁▁▁
Sale_Type_ConLI 0 1 0 1 -0.05 -0.05 -0.05 -0.05 18.27 ▇▁▁▁▁
Sale_Type_ConLw 0 1 0 1 -0.05 -0.05 -0.05 -0.05 18.27 ▇▁▁▁▁
Sale_Type_CWD 0 1 0 1 -0.07 -0.07 -0.07 -0.07 13.94 ▇▁▁▁▁
Sale_Type_New 0 1 0 1 -0.30 -0.30 -0.30 -0.30 3.35 ▇▁▁▁▁
Sale_Type_Oth 0 1 0 1 -0.05 -0.05 -0.05 -0.05 19.74 ▇▁▁▁▁
Sale_Type_VWD 0 1 0 1 -0.02 -0.02 -0.02 -0.02 48.40 ▇▁▁▁▁
Sale_Type_WD. 0 1 0 1 -2.51 0.40 0.40 0.40 0.40 ▁▁▁▁▇
Sale_Condition_AdjLand 0 1 0 1 -0.07 -0.07 -0.07 -0.07 13.94 ▇▁▁▁▁
Sale_Condition_Alloca 0 1 0 1 -0.09 -0.09 -0.09 -0.09 11.70 ▇▁▁▁▁
Sale_Condition_Family 0 1 0 1 -0.13 -0.13 -0.13 -0.13 7.69 ▇▁▁▁▁
Sale_Condition_Normal 0 1 0 1 -2.13 0.47 0.47 0.47 0.47 ▂▁▁▁▇
Sale_Condition_Partial 0 1 0 1 -0.30 -0.30 -0.30 -0.30 3.30 ▇▁▁▁▁

Aplica-se então cross validation na base de treino. Além disso, serão criadas amostras bootstrap para ajuste dos hiperparâmetros de alguns dos modelos utilizados.

set.seed(123)
(cv_splits <- vfold_cv(train_baked, v = 5, strata = Sale_Price))
## #  5-fold cross-validation using stratification 
## # A tibble: 5 x 2
##   splits             id   
##   <named list>       <chr>
## 1 <split [1.9K/470]> Fold1
## 2 <split [1.9K/470]> Fold2
## 3 <split [1.9K/469]> Fold3
## 4 <split [1.9K/468]> Fold4
## 5 <split [1.9K/468]> Fold5
set.seed(123)
(ames_boot <- bootstraps(train_baked, times = 10, strata = Sale_Price))
## # Bootstrap sampling using stratification 
## # A tibble: 10 x 2
##    splits             id         
##    <named list>       <chr>      
##  1 <split [2.3K/865]> Bootstrap01
##  2 <split [2.3K/857]> Bootstrap02
##  3 <split [2.3K/890]> Bootstrap03
##  4 <split [2.3K/858]> Bootstrap04
##  5 <split [2.3K/872]> Bootstrap05
##  6 <split [2.3K/864]> Bootstrap06
##  7 <split [2.3K/855]> Bootstrap07
##  8 <split [2.3K/856]> Bootstrap08
##  9 <split [2.3K/844]> Bootstrap09
## 10 <split [2.3K/876]> Bootstrap10

4. Modelos

Antes de ajustar os modelos, é necessário especificar o pacote usado, modo e parâmetros. Depois disso, é feita a modelagem nos diferentes folds criados por cross validation para avaliar seu desempenho.

4.1. Regressão Linear

4.1.1. Sem seleção stepwise

lm_spec <- linear_reg() %>% 
  set_engine('lm')

#fit nos cv folds
lm_res <- fit_resamples(Sale_Price ~ .,
               lm_spec,
               cv_splits,
               control = control_resamples(save_pred = TRUE))

Resumo do modelo ajustado:

(lm_res %>% 
  collect_metrics())
## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.335     5  0.0336
## 2 rsq     standard   0.887     5  0.0197

Modelo ajustado na base completa, resumo e importância das variáveis:

lm_fit <- lm_spec %>% 
      fit(Sale_Price ~.,
      data = train_baked)

lm_fit %>% 
  summary()
##         Length Class      Mode   
## lvl      0     -none-     NULL   
## spec     5     linear_reg list   
## fit     12     lm         list   
## preproc  1     -none-     list   
## elapsed  5     proc_time  numeric
lm_fit %>% 
  tidy()
## # A tibble: 187 x 5
##    term          estimate std.error statistic  p.value
##    <chr>            <dbl>     <dbl>     <dbl>    <dbl>
##  1 (Intercept)  -1.28e-14   0.00580 -2.21e-12 1.00e+ 0
##  2 Lot_Frontage  1.95e- 2   0.00720  2.71e+ 0 6.74e- 3
##  3 Lot_Area      4.44e- 2   0.00776  5.72e+ 0 1.20e- 8
##  4 Lot_Shape    -5.80e- 3   0.00738 -7.86e- 1 4.32e- 1
##  5 Land_Contour -8.18e- 3   0.00803 -1.02e+ 0 3.08e- 1
##  6 Utilities     1.75e- 3   0.00712  2.47e- 1 8.05e- 1
##  7 Land_Slope    2.35e- 3   0.00829  2.84e- 1 7.77e- 1
##  8 Overall_Qual  1.48e- 1   0.0128   1.15e+ 1 7.32e-30
##  9 Overall_Cond  7.26e- 2   0.00852  8.52e+ 0 3.07e-17
## 10 Year_Built    1.54e- 1   0.0221   6.98e+ 0 4.05e-12
## # … with 177 more rows

4.1.2 Com seleção stepwise backward

O pacote tidymodels não possui uma integração pronta com o pacote leaps para realizar stepwise selection. Devido a isso, o procedimento será feito, mas o modelo não será considerado na avaliação final. É possível criar essa integração pelo pacote parsnip.

library(leaps)

stepb <- regsubsets(Sale_Price ~ ., data = train_baked, nvmax = 10, 
                          method = "backward")
## Warning in leaps.setup(x, y, wt = wt, nbest = nbest, nvmax = nvmax, force.in =
## force.in, : 3 linear dependencies found
## Reordering variables and trying again:
  resumo_stepb <- summary(stepb)
  #resumo_stepb$outmat
  
  coef(stepb, id = which.min(resumo_stepb$bic))
##       (Intercept)      Overall_Qual        Exter_Qual       Bsmt_Unf_SF 
##     -2.632578e-17      3.166017e-01      2.124920e-01     -1.651529e-01 
##     Total_Bsmt_SF      First_Flr_SF     Second_Flr_SF Misc_Feature_Shed 
##      2.491390e-01      3.179770e-01      3.014907e-01     -3.175356e-03 
## Misc_Feature_TenC     Sale_Type_Con   Sale_Type_ConLD     Sale_Type_WD. 
##     -9.898401e-03      3.379658e-03     -1.065269e-02     -5.593790e-02
  plot(stepb, scale = "adjr2")

4.2. Lasso

Nesse modelo, é utilizada a função tune() para encontrar o melhor lambda (penalty) nas amostras bootstrap.

lasso_spec <- linear_reg(penalty = tune(), mixture = 1) %>%
  set_engine("glmnet", standardize = FALSE) #false, pois a base já foi pré-processada anteriormente 

lambda_grid <- grid_regular(penalty(), levels = 500)

doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(123)
lasso_grid <- tune_grid(Sale_Price ~.,
  model = lasso_spec,
  resamples = ames_boot,
  grid = lambda_grid)

Avaliação dos parâmetros - lambda para o menor erro quadrado médio (rmse):

lasso_grid %>%
  collect_metrics()
## # A tibble: 1,000 x 6
##     penalty .metric .estimator  mean     n std_err
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
##  1 1.00e-10 rmse    standard   0.370    10 0.0128 
##  2 1.00e-10 rsq     standard   0.870    10 0.00796
##  3 1.05e-10 rmse    standard   0.370    10 0.0128 
##  4 1.05e-10 rsq     standard   0.870    10 0.00796
##  5 1.10e-10 rmse    standard   0.370    10 0.0128 
##  6 1.10e-10 rsq     standard   0.870    10 0.00796
##  7 1.15e-10 rmse    standard   0.370    10 0.0128 
##  8 1.15e-10 rsq     standard   0.870    10 0.00796
##  9 1.20e-10 rmse    standard   0.370    10 0.0128 
## 10 1.20e-10 rsq     standard   0.870    10 0.00796
## # … with 990 more rows
lasso_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
    ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  labs(x = 'Log(lambda)') +
  theme(legend.position = "none")
## Warning: Removed 5 rows containing missing values (geom_errorbar).
## Warning: Removed 4 rows containing missing values (geom_path).

(lasso_lowest_rmse <- lasso_grid %>%
  select_best("rmse", maximize = FALSE))
## # A tibble: 1 x 1
##   penalty
##     <dbl>
## 1 0.00359

Modelo Lasso final: O pacote tidymodels permite finalizar o modelo com o melhor parâmetro encontrado. Após isso, o modelo é aplicado nas amostras criadas com cross-validatione também ajustada na base de treino completa.

lasso_final <- lasso_spec %>% 
  finalize_model(parameters = lasso_lowest_rmse)

lasso_res <- fit_resamples(Sale_Price ~ .,
               lasso_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))
    
(lasso_res %>% 
  collect_metrics())
## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.340     5  0.0289
## 2 rsq     standard   0.883     5  0.0171
lasso_fit <- lasso_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

Infelizmente essa forma de ajustar os modelos não permite a visualização das variáveis pela biblioteca plotmo, pois a classe criada pelo tidymodels não é reconhecida. É possível, no entanto, montar o gráfico, mas nesse caso específico, devido à quantidade de variáveis, a visualização fica prejudicada.

#library(plotmo)
#plot(lasso_fit, xvar = 'lambda', label = FALSE)

#?plot_glmnet

### Fazer fit na base inteira
#plot_glmnet(lasso_fit)

lasso_fit %>% 
  tidy %>% 
  mutate(log.lambda = log(lambda)) %>% 
  #filter(step <=100) %>% 
  #filter(estimate < 0.2) %>% 
  ggplot(aes(log.lambda, estimate, color = term)) +
  geom_line(show.legend = FALSE)

Será utilizado, portanto, o pacote vip, que possui funcionalidades semelhantes e permite identificar as variáveis mais importantes no modelo.

#plot(lasso.coef[lasso.coef != 0])
lasso_var <- lasso_fit %>%
  vi(lambda = lasso_lowest_rmse$penalty) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct))

#Verificaçãoda seleção das variáveis
lasso_var %>% 
  count(Importance_pct != 0)
## # A tibble: 2 x 2
##   `Importance_pct != 0`     n
##   <lgl>                 <int>
## 1 FALSE                     4
## 2 TRUE                    185

Aparentemente o modelo removeu apenas 4 variáveis, o que não ajuda muito nesse caso.

#variáveis que mais impactam no preço  
lasso_var %>% 
  filter(Importance_pct > 0.05) %>% 
  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
  geom_col()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

#variáveis que impactam negativamente no preço
#lasso_var %>% 
#  filter(Importance < 0) %>% 
#  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
#  geom_col()+
#  scale_y_continuous(labels = scales::percent_format())+
#  coord_flip()

Observa-se que as principais variáveis são da mesma categoria e não fazem muito sentido se compararmos com o que foi visto na análise exploratória. Talvez rodar o modelo novamente com outros parâmetros e pré-processamento deva ser considerado. Foi filtrado o percntual acima de 5% para o gráfico não ficar muito poluído.

4.3. Ridge Regression

Assim, como no Lasso, o lambda será otimizado.

ridge_spec <- linear_reg(penalty = tune(), mixture = 0) %>% #mixture = 0 para ridge regression
  set_engine("glmnet", standardize = FALSE)

ridge_lambda_grid <- grid_regular(penalty(), levels = 1000)

doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(2020)
ridge_grid <- tune_grid(Sale_Price ~.,
  model = ridge_spec,
  resamples = ames_boot,
  grid = ridge_lambda_grid
)
ridge_grid %>%
  collect_metrics()
## # A tibble: 2,000 x 6
##     penalty .metric .estimator  mean     n std_err
##       <dbl> <chr>   <chr>      <dbl> <int>   <dbl>
##  1 1.00e-10 rmse    standard   0.355    10 0.0120 
##  2 1.00e-10 rsq     standard   0.879    10 0.00728
##  3 1.02e-10 rmse    standard   0.355    10 0.0120 
##  4 1.02e-10 rsq     standard   0.879    10 0.00728
##  5 1.05e-10 rmse    standard   0.355    10 0.0120 
##  6 1.05e-10 rsq     standard   0.879    10 0.00728
##  7 1.07e-10 rmse    standard   0.355    10 0.0120 
##  8 1.07e-10 rsq     standard   0.879    10 0.00728
##  9 1.10e-10 rmse    standard   0.355    10 0.0120 
## 10 1.10e-10 rsq     standard   0.879    10 0.00728
## # … with 1,990 more rows
ridge_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
  ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  theme(legend.position = "none")

(ridge_lowest_rmse <- ridge_grid %>%
  select_best("rmse", maximize = FALSE))
## # A tibble: 1 x 1
##   penalty
##     <dbl>
## 1   0.117

Modelo Final:

ridge_final <- ridge_spec %>% 
  finalize_model(parameters = ridge_lowest_rmse)

ridge_res <- fit_resamples(Sale_Price ~ .,
               ridge_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))

(ridge_res %>% 
  collect_metrics())
## # A tibble: 2 x 5
##   .metric .estimator  mean     n std_err
##   <chr>   <chr>      <dbl> <int>   <dbl>
## 1 rmse    standard   0.340     5  0.0273
## 2 rsq     standard   0.884     5  0.0159
ridge_fit <- ridge_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

ridge_fit %>% 
  summary
##         Length Class      Mode   
## lvl      0     -none-     NULL   
## spec     6     linear_reg list   
## fit     12     elnet      list   
## preproc  5     -none-     list   
## elapsed  5     proc_time  numeric
ridge_fit %>% 
  tidy %>% 
  mutate(log.lambda = log(lambda)) %>% 
  #filter(step <=100) %>% 
  #filter(estimate < 0.2) %>% 
  ggplot(aes(log.lambda, estimate, color = term)) +
  geom_line(show.legend = FALSE)

#plot(lasso.coef[lasso.coef != 0])
var_ridge <- ridge_fit %>%
  vi(lambda = ridge_lowest_rmse$penalty) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct))

#variáveis que mais impactam no preço  
var_ridge %>% 
  filter(Importance_pct > 0.20) %>% 
  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
  geom_col()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

#variáveis que impactam negativamente no preço
#var_ridge %>% 
#  filter(Importance < 0) %>% 
#  ggplot(aes(Variable,Importance_pct, fill = Sign)) +
#  geom_col()+
#  scale_y_continuous(labels = scales::percent_format())+
#  coord_flip()

Apesar de não zerar o coeficiente de nenhuma variável, o modelo ridge selecionou varíveis que estão mais correlacionadas com a variável resposta. Além disso, observa-se o impacto negativo de algumas das principais selecionadas (mais de 50%).

4.4. Bagging

Assim como a seleção stepwise, também não há uma função específica no tidymodels para bagging. Seria possível através do crescimento de várias árvores, mas para simplificar, será utilizado o pacote ipred que realiza esse procedimento muito mais rápido.

library(ipred)

n_trees = seq(50,2200,500)

#função para calcular o eqm
f_eqm <- function(model, y = train_baked$Sale_Price){
          sum((y - round(model$y))^2) / length(model$y)
        }

bag_res <- tibble(trees = n_trees,
                  oob = NA,
                  rmse = NA)

for (i in 1:nrow(bag_res)){
  
  #print(i)
  
  set.seed(123)
  bag_fit <- ipred::bagging(Sale_Price ~ ., data = train_baked, coob = TRUE,  nbagg = n_trees[i])
  bag_res$oob[i] <- bag_fit$err
  bag_res$rmse[i] <- f_eqm(bag_fit)
}

A partir dos resultados computados, avalia-se o erro out-of-bag para diferentes valores de árvores:

bag_res %>% 
  ggplot(aes(trees,oob))+
  geom_line()

Como podemos observar, mesmo sendo uma diferença pequena, o erro estabiliza a partir de 1000 “bagged” árvores. Com isso, temos o modelo final.

bag_fit <- ipred::bagging(Sale_Price ~ ., data = train_baked, coob = TRUE, nbagg = 500)

4.5. Random Forest

Aqui será feita a otimização do parâmetro mtry e o número de árvores. Para otimizar o processamento, iniciaremos com um grid apenas para o mtry, para depois testar o número de árvores.

p <- ncol(train_baked) - 1 #total de variáveis preditoras

rf_tune_spec <- rand_forest(mode = "regression",
                            mtry = tune(), #p/3 = 189
                            trees = 500) %>% 
                set_engine("ranger", importance = "permutation")


(rf_grid <- grid_regular(
  mtry(range = c(p/5, 100)),
  #trees(range = c(500, 2000)),
  levels = 5))
## # A tibble: 5 x 1
##    mtry
##   <int>
## 1    37
## 2    53
## 3    68
## 4    84
## 5   100
doParallel::registerDoParallel() #processamento em paralelo para otimizar

set.seed(2020)
rf_tune <- tune_grid(Sale_Price ~.,
                      model = rf_tune_spec,
                      resamples = ames_boot,
                      grid = rf_grid)

Avaliando os resultados:

rf_tune %>% 
  collect_metrics()
## # A tibble: 10 x 6
##     mtry .metric .estimator  mean     n std_err
##    <int> <chr>   <chr>      <dbl> <int>   <dbl>
##  1    37 rmse    standard   0.326    10 0.0117 
##  2    37 rsq     standard   0.904    10 0.00630
##  3    53 rmse    standard   0.325    10 0.0114 
##  4    53 rsq     standard   0.903    10 0.00619
##  5    68 rmse    standard   0.325    10 0.0115 
##  6    68 rsq     standard   0.903    10 0.00632
##  7    84 rmse    standard   0.326    10 0.0115 
##  8    84 rsq     standard   0.901    10 0.00642
##  9   100 rmse    standard   0.327    10 0.0117 
## 10   100 rsq     standard   0.900    10 0.00662
rf_tune%>%
  collect_metrics() %>%
  select(mean, mtry, .metric) %>%
  filter(.metric == 'rmse') %>% 
  #pivot_longer(min_n:mtry,
  #  values_to = "value",
  #  names_to = "parameter"
  #) %>%
  ggplot(aes(mtry, mean, color = .metric)) +
  geom_point(show.legend = TRUE) #+

  #facet_wrap(~parameter, scales = "free_x") +
  #labs(x = NULL, y = "Value")

rf_lowest_rmse <- rf_tune %>%
  select_best("rmse", maximize = FALSE)

best_mtry = rf_lowest_rmse$mtry

Observa-se que mtry ótimo está perto de 50. Acima disso, o modelo provavelmente está overfitting. Com o melhor mtry, vamos avaliar o erro considerando o número de árvores e finalizar o modelo com os melhores parâmetros.

  • Melhor número de árvores
rf_spec <- rand_forest(mode = "regression",
                            mtry = best_mtry,
                            trees = tune()) %>% 
                set_engine("ranger", importance = "permutation")


rf_grid <- grid_regular(
  #mtry(range = c(p/4,100)),
  trees(range = c(500, 2000)),
  levels = 4)

doParallel::registerDoParallel()

set.seed(2020)
rf_tune <- tune_grid(Sale_Price ~.,
                      model = rf_spec,
                      resamples = ames_boot,
                      grid = rf_grid)

rf_tune %>% 
  collect_metrics()
## # A tibble: 8 x 6
##   trees .metric .estimator  mean     n std_err
##   <int> <chr>   <chr>      <dbl> <int>   <dbl>
## 1   500 rmse    standard   0.325    10 0.0115 
## 2   500 rsq     standard   0.903    10 0.00629
## 3  1000 rmse    standard   0.324    10 0.0114 
## 4  1000 rsq     standard   0.904    10 0.00621
## 5  1500 rmse    standard   0.324    10 0.0113 
## 6  1500 rsq     standard   0.904    10 0.00624
## 7  2000 rmse    standard   0.324    10 0.0114 
## 8  2000 rsq     standard   0.904    10 0.00625
rf_tune %>%
  collect_metrics() %>%
  select(mean, trees, .metric) %>%
  filter(.metric == 'rmse') %>% 
  #pivot_longer(min_n:mtry,
  #  values_to = "value",
  #  names_to = "parameter") %>%
  ggplot(aes(trees, mean, color = .metric)) +
  geom_point(show.legend = TRUE)

  rf_lowest_rmse <- rf_tune %>%
    select_best("rmse", maximize = FALSE)
  • Modelo Final:
rf_final <-  rf_spec %>%
    finalize_model(parameters = rf_lowest_rmse)


rf_res <- fit_resamples(Sale_Price ~ .,
               rf_final,
               cv_splits,
               control = control_resamples(save_pred = TRUE))


rf_fit <- rf_final %>% 
      fit(Sale_Price ~.,
      data = train_baked)

rf_fit %>% 
  summary
##         Length Class       Mode   
## lvl      0     -none-      NULL   
## spec     6     rand_forest list   
## fit     15     ranger      list   
## preproc  1     -none-      list   
## elapsed  5     proc_time   numeric

Importância das variáveis

vi(rf_fit) %>% 
  mutate(Importance_pct = abs(Importance)/max(abs(Importance))) %>% 
  mutate(Variable = fct_reorder(Variable, Importance_pct)) %>% 
  filter(Importance_pct > 0.05) %>% 
  ggplot(aes(Variable, Importance_pct)) +
  geom_point()+
  scale_y_continuous(labels = scales::percent_format())+
  coord_flip()

5. Comparação dos Modelos

No gráfico abaixo, é possível observar o erro quadrado médio (rmse) e o R-quadrado (rsq) para cada um dos modelos, exceto para bagging, que foi ajustado diretamente na base completa.

lm_res %>% 
  select(id, .metrics) %>% 
  unnest(.metrics) %>% 
  mutate(model = "linear regression") %>% 
  bind_rows(lasso_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "lasso")) %>% 
   bind_rows(ridge_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "ridge")) %>% 
  #bind_rows(bag_res %>% 
   #         select(id, .metrics) %>% 
    #        unnest(.metrics) %>% 
     #       mutate(model = "bagging")) %>% 
  bind_rows(rf_res %>% 
            select(id, .metrics) %>% 
            unnest(.metrics) %>% 
            mutate(model = "random forest")) %>% 
  ggplot(aes(id, .estimate, group = model, color = model)) + 
    geom_point(size = 1.5) + 
    facet_wrap(~.metric) + 
    coord_flip()

Para ter uma ideia melhor, segue a comparação das métricas nas bases de treino e teste completas.

results_train <- lm_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'lm') %>% 
  bind_rows(lasso_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'lasso')) %>%
  bind_rows(ridge_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'ridge')) %>%
  bind_rows(tibble(.pred = predict(bag_fit, newdata = train_baked),
                    truth = train_baked$Sale_Price,
                    model = 'bagging')) %>% 
  bind_rows(rf_fit %>% 
  predict(new_data = train_baked) %>% 
  mutate(truth = train_baked$Sale_Price,
         model = 'random forest'))

results_train %>% 
  group_by(model) %>% 
  rmse(truth = truth, estimate = .pred)
## # A tibble: 5 x 4
##   model         .metric .estimator .estimate
##   <chr>         <chr>   <chr>          <dbl>
## 1 bagging       rmse    standard       0.393
## 2 lasso         rmse    standard       0.290
## 3 lm            rmse    standard       0.270
## 4 random forest rmse    standard       0.128
## 5 ridge         rmse    standard       0.295

Pela base de treino, observa-se que o melhor modelo foi random forest, seguido pela regressão linear. Para confirmar, vamos observar o comportamento na base teste:

results_test <- lm_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'lm') %>% 
  bind_rows(lasso_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'lasso')) %>%
  bind_rows(ridge_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'ridge')) %>%
  bind_rows(tibble(.pred = predict(bag_fit, newdata = test_baked),
                    truth = test_baked$Sale_Price,
                    model = 'bagging')) %>% 
  bind_rows(rf_fit %>% 
  predict(new_data = test_baked) %>% 
  mutate(truth = test_baked$Sale_Price,
         model = 'random forest'))
## Warning in predict.lm(object = object$fit, newdata = new_data, type =
## "response"): prediction from a rank-deficient fit may be misleading
results_test %>% 
  group_by(model) %>% 
  rmse(truth = truth, estimate = .pred)
## # A tibble: 5 x 4
##   model         .metric .estimator .estimate
##   <chr>         <chr>   <chr>          <dbl>
## 1 bagging       rmse    standard       0.481
## 2 lasso         rmse    standard       0.503
## 3 lm            rmse    standard       0.557
## 4 random forest rmse    standard       0.375
## 5 ridge         rmse    standard       0.479

Na base de teste, tivemos o mesmo resultado: o modelo que melhor performou foi random forest. No entanto, o segundo melhor modelo foi bagging, que teve a pior performance na base de treino. Ridge e Lasso também tiveram a performance invertida. A regressão linear, no entanto, foi o modelo que performou pior nesse caso.

6. Conclusão

A floresta aleatória foi o modelo que melhor performou, tanto na base treino, como na base teste e, portanto, o modelo selecionado como melhor. Conforme avaliado na análise exploratória, as principais variáveis preditoras estavam altammente correlacionadas, mas não necessariamente no topo da lista. As 3 principais foram: Overall_Quality, Gr_Liv_Area e Year_Built.